{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![cant.png](figures/helmholtz.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Method of Manufactured solutions and solve the Helmholtz equation. \n",
    "\n",
    "Problem taken from http://www.dealii.org/8.4.1/doxygen/deal.II/step_7.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuAFEM\n",
    "using KrylovMethods\n",
    "using Tensors\n",
    "const ∇ = Tensors.gradient\n",
    "const Δ = Tensors.hessian;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = generate_grid(Quadrilateral, (150, 150))\n",
    "addnodeset!(grid, \"dirichlet_boundary\", x -> x[1] ≈ 1 ||  x[2] ≈ 1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 2\n",
    "ip = Lagrange{dim, RefCube, 1}()\n",
    "qr = QuadratureRule{dim, RefCube}(2)\n",
    "qr_face = QuadratureRule{dim-1, RefCube}(2)\n",
    "cellvalues = CellScalarValues(qr, ip);\n",
    "facevalues = FaceScalarValues(qr_face, ip);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DofHandler\n",
       "  Fields:\n",
       "    u dim: 1\n",
       "  Total dofs: 22801\n",
       "  Dofs per cell: 4"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dh = DofHandler(grid)\n",
    "push!(dh, :u, 1)\n",
    "close!(dh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function u_ana{T}(x::Vec{2, T}) \n",
    "    xs = (Vec{2}((-0.5,  0.5)),\n",
    "          Vec{2}((-0.5, -0.5)),\n",
    "          Vec{2}(( 0.5,  -0.5)))\n",
    "    σ = 1/8\n",
    "    s = zero(eltype(x))\n",
    "    for i in 1:3\n",
    "        s += exp(- norm(x - xs[i])^2 / σ^2)\n",
    "    end\n",
    "    return max(1e-15 * one(T), s) # Denormals, be gone\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dbc = DirichletBoundaryConditions(dh)\n",
    "add!(dbc, :u, getnodeset(grid, \"dirichlet_boundary\"), (x,t) -> u_ana(x))\n",
    "close!(dbc)\n",
    "update!(dbc, 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = create_sparsity_pattern(dh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function doassemble{dim}(cellvalues::CellScalarValues{dim}, facevalues::FaceScalarValues{dim},\n",
    "                         K::SparseMatrixCSC, dh::DofHandler)\n",
    "    b = 1.0\n",
    "    f = zeros(ndofs(dh))\n",
    "    assembler = start_assemble(K, f)\n",
    "    \n",
    "    n_basefuncs = getnbasefunctions(cellvalues)\n",
    "    global_dofs = zeros(Int, ndofs_per_cell(dh))\n",
    "\n",
    "    fe = zeros(n_basefuncs) # Local force vector\n",
    "    Ke = zeros(n_basefuncs, n_basefuncs) # Local stiffness mastrix\n",
    "\n",
    "    @inbounds for (cellcount, cell) in enumerate(CellIterator(dh))\n",
    "        fill!(Ke, 0)\n",
    "        fill!(fe, 0)\n",
    "        coords = getcoordinates(cell)\n",
    "\n",
    "        reinit!(cellvalues, cell)\n",
    "        for q_point in 1:getnquadpoints(cellvalues)\n",
    "            dΩ = getdetJdV(cellvalues, q_point)\n",
    "            \n",
    "            coords_qp = spatial_coordinate(cellvalues, q_point, coords)\n",
    "            f_true = trace(hessian(u_ana, coords_qp)) + u_ana(coords_qp)\n",
    "            for i in 1:n_basefuncs\n",
    "                δu = shape_value(cellvalues, q_point, i)\n",
    "                ∇δu = shape_gradient(cellvalues, q_point, i)\n",
    "                fe[i] += (δu * f_true) * dΩ\n",
    "                for j in 1:n_basefuncs\n",
    "                    u = shape_value(cellvalues, q_point, j)\n",
    "                    ∇u = shape_gradient(cellvalues, q_point, j)\n",
    "                    Ke[i, j] += (∇δu ⋅ ∇u + δu * u) * dΩ\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "        \n",
    "        for face in 1:nfaces(cell)\n",
    "            if onboundary(cell, face) && \n",
    "                   ((cellcount, face) ∈ getfaceset(grid, \"left\") || \n",
    "                    (cellcount, face) ∈ getfaceset(grid, \"bottom\"))\n",
    "                reinit!(facevalues, cell, face)\n",
    "                for q_point in 1:getnquadpoints(facevalues)\n",
    "                    coords_qp = spatial_coordinate(facevalues, q_point, coords)\n",
    "                    n = getnormal(facevalues, q_point)\n",
    "                    g = gradient(u_ana, coords_qp) ⋅ n\n",
    "                    dΓ = getdetJdV(facevalues, q_point)\n",
    "                    for i in 1:n_basefuncs\n",
    "                        δu = shape_value(facevalues, q_point, i)\n",
    "                        fe[i] += (δu * g) * dΓ\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "   \n",
    "        celldofs!(global_dofs, cell)\n",
    "        assemble!(assembler, global_dofs, fe, Ke)\n",
    "    end\n",
    "    return K, f\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K, f = doassemble(cellvalues, facevalues, K, dh);\n",
    "apply!(K, f, dbc)\n",
    "u = cholesky(Symmetric(K)) \\ f;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{String,1}:\n",
       " \"helmholtz.vtu\""
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vtkfile = vtk_grid(\"helmholtz\", dh)\n",
    "vtk_point_data(vtkfile, dh, u)\n",
    "vtk_save(vtkfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Helmholtz successful\n"
     ]
    }
   ],
   "source": [
    "Base.Test.@test maximum(u) ≈ 0.05637592090022005\n",
    "println(\"Helmholtz successful\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "julia-0.6",
   "display_name": "Julia 0.6.0",
   "language": "julia"
  },
  "language_info": {
   "mimetype": "application/julia",
   "file_extension": ".jl",
   "version": "0.6.0",
   "name": "julia"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
